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The prediction of the quasi- static response of industrial 
laminate structures requires to use fine descriptions of the ma- 
terial, especially when debonding is involved. Even when 
modeled at the mesoscale, the computation of these struc- 
tures results in very large numerical problems. In this pa- 
per, the exact mesoscale solution is sought using parallel it- 
erative solvers. The LaTIn-based mixed domain decompo- 
sition method makes it very easy to handle the complex de- 
scription of the structure; moreover the provided multiscale 
features enable us to deal with numerical difficulties at their 
natural scale; we present the various enhancements we devel- 
oped to ensure the scalability of the method. An extension of 
the method designed to handle instabilities is also presented. 

1 Introduction 

Since the early 1980's, a very large number of studies has 
been conducted on the prediction of debonding in compos- 
ite laminates, resulting in better understanding of the failure 
processes of composites. As a result, micromechanics-based 
models have been shown to enable accurate predictions of the 
debonding phenomenon. 

The industrialists' wish to replace expensive experiments 
by virtual tests for the design of their composite structures 
has brought about a new issue. Currently, the analysis of in- 
dustrial problems using previously referenced micromodels is 
infeasible because the memory size and computing time re- 
quirements are far too large. Therefore, most applications to 
predict the initiation and propagation of debonding in lam- 
inates rely on mesoscale modeling 17] [TSl which are also 
rooted in the analysis of micromechanics phenomena. In this 
paper we retain the model described in IS |3l ; it describes the 
behavior of the laminates distinctly in the plies, 3D entities, 
and in the interfaces, 2D entities, the debonding ability be- 
ing localized in the interfaces and handled through a cohesive 
behavior. 



Even for small test cases, the numerical problem resulting 
from the mesomodeling of laminate structures is huge. Nev- 
ertheless, the latest advances in domain decomposition and 
multiscale methods provide powerful tools enabling the cal- 
culation of laminate structures of reasonable size. We can 
distinguish two families of solvers able to handle such large 
numerical problems. The first one consists in using a non- 
linear homogenized strategy |[29l [TTl [Tol coupled with local 
enrichment methods ['^[ElllSlEllinil^gl. In this paper we 
focus on the prediction of debonding in process zones within 
the structure of the laminate, in which we consider that the 
debonding processes leading to final failure can be circum- 
scribed. Though, within these potentially large process zones, 
no low gradient zone can be identified a priori. Thus, the 
enrichment-based strategy might be difficult to use. Instead, 
we wish to use the second family of solvers for large numeri- 
cal problems, which consists in computing the exact mesoso- 
lution everywhere in the process zones, using parallel iterative 
solvers Li,22llI21[Il[l2l[l. Coupling the solution with a re- 
duced model such as a plate model |[3T||4l will be the subject 
of further work, and is not dealt with in this paper. 

The mixed domain decomposition strategy described in 
1 21 1 uses an original concept which consists in splitting the 
structure in volume substructures separated by surface inter- 
faces. Consequently, the reference problem resulting from the 
chosen mesomodeling is very naturally substructured, the co- 
hesive interfaces of the model being handled within the inter- 
faces of the domain decomposition method. This idea is de- 
veloped in Section ([2]). Furthermore, the resolution of the sub- 
structured problem by a LaTIn iterative solver exhibits very 
interesting numerical properties: the nonlinearities are dealt 
with through local problems and very few re-assembling steps 
are required. The incremental micro-macro LaTIn algorithm 
is presented in Section \22) without any improvement. As 
shown in Section \23) several numerical difficulties are en- 
countered when directly applying the method. The core of the 
paper (Sections [3] to [6]) presents and assesses improvements 
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and adaptations of the method to efficiently handle delamina- 
tion computations. 

The substructured LaTIn method is parametrized by two in- 
terface search direction operators. Optimal values have been 
determined and practical values have been assessed for many 
mechanical problems (perfect interfaces, contact with or with- 
out friction), though for cohesive interfaces an adaptation is 
required in order the method to be efficient (and in some case 
feasible), as explained in Section ([3]). 

The method is granted a multiscale nature by the separa- 
tion of a macroscopic part and a microscopic part of the in- 
terface fields. This separation, coupled with continuity condi- 
tions, leads to the construction of an automatic homogeniza- 
tion procedure. This concept has been successfully applied 
in the cases of crack propagation problems in homogeneous 
media in |15|, of composite structures modeled at the mi- 
croscale 1 19], and of multiscale problems in time and space 
II20I . Though, in the case of singularities resulting from the 
crack tips being localized on the interfaces of the domain de- 
composition strategy, the homogeneous solution is too poor 
to represent accurately the solution. It results in a loss of ex- 
tensibility of the strategy. Therefore, we enhance the method 
with a specific technique for the calculation of the quantities 
in the process zone with increased accuracy, which results in 
a significant improvement of the convergence rate, which is 
the topic of Section (|4]). 

A consequence to the sub structuring naturally introduced 
to solve the reference problem is that the number of substruc- 
tures required to solve large delamination problems becomes 
huge. Hence, the macroscopic problem can not be solved us- 
ing direct solvers. The introduction of a third-level problem is 
required to quickly propagate the very-high-wavelength part 
of the solution. This is achieved by solving the second-scale 
problem using the balancing domain decomposition method 
described in |22|, as explained in Section ([5]). 

When trying to predict the very final residual strength of 
the structure, which is of the industrialist's interest, one has 
to deal with instabilities and limit-points problems resulting 
from the local softening behavior. Hence, an adaptation of 
the three-scale domain decomposition method to arc-length- 
type algorithms with local control of the loading amplitude 
has been developed, which is the subject of Section ([6]). 

2 The reference problem and strategy 

2.1 Reference problem and substructuring 

The laminate structure E occupying the domain is made 
out of Np adjacent plies P occupying Domain Q.p, separated 
by Np — I cohesive interfaces Ippf (see Figure ([T])). An ex- 
ternal traction field (respectively a displacement field U^) 




Figure 1 : Reference problem 



is applied to the structure on Part dQf (respectively dQu) of 
the boundary dQ, of Domain Q, = [jpQ,p. The normal to the 
boundary dQ,p of Ply P, external to P, is Hp. The volume 
force is denoted f^. Let Up be the displacement field, the 
Cauchy stress tensor and the symmetric part of the dis- 
placement gradient in Ply P. 

The simulation is performed under the assumption of small 
perturbations and the evolution over time is considered to be 
quasi- static and isothermal. The problem is solved using an 
implicit time integration scheme. At each time step of the 
analysis, the reference equilibrium problem is: 

Find Sref = {sp)p^^, where sp = {o^^Up), which verifies 
the following equations: 

• Kinematic admissibility on dClp D dQ^u- 

Up = U^ (1) 

• Global equilibrium of Structure E: 

V(m/)/>6e e '^i'^ X ... X 

pJap-^ pJdnndnf 

+ 11 /, , ^pnpMp*dr = 

P pT^pJdapndap,—^ — ^ 

(2) 

• Linear orthotropic behavior of the plies: 

at each point of Q,p^ —p~ ^^(^p) 

• Constitutive equation of the interfaces: 

at each point of Ipp/ , ^pp/ ( [w]^ , GpUp ) = (4) 

The gap of displacement [u]^ of Interface Ippf such that 
P <P' has arbitrarily been introduced as [u]^ = Upi — Up. The 
operator s^ppt establishes a relation between the primal in- 
terface unknown [w]^, and the dual interface unknown a^^p, 
which reads : 

gp'Hp = Kppf .[u[p (5) 
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The expression of the local stiffness operator K^^, of Interface 
Ippi can be made explicit in the basis {np-,t_i^t_2-,) (see Figure 
0): 

( (l-h^{[iAp-np)d,)kl 

Epp'= (i-^i)^? 

V {l-d2)k^ 

is here the positive indicator function. 






Up 






A_^^ 




Figure 2: The mesomodel entities 

The local damage variables di are introduced into the inter- 
face model in order to simulate its softening behavior when 
the structure is loaded. Their values range from (healthy in- 
terface point) to 1 (completely damaged interface point). The 
parameters di are related to the local energy release rates Yi 
of the interface's degradation modes. Denoting the surface 
strain energy of the cohesive interface, 

(6) 

grf is here the surface strain energy of the cohesive interface. 

We assume that the damage variables are functions of a sin- 
gle quantity: the maximum Y^^ over time of a combination of 
the energy release rates i^i^, x <t: 



Yi 



ddi 



where 



a 

3|T- 



Y\t = max(^<,) 

Thus, the evolution laws are: 

d\= d2= d^ = w(Y) 

where in general w(Y) = 



(7) 



(8) 
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This model has the advantage of using a single damage 
variable to handle several macroscopic delamination modes 
of the interface (traction along Hp and shear along t_i and 
t_2)- However, when setting Parameters Y\ and 72 to identified 
physical values such that 7i ^ 72 / 1, the energy dissipated 
due to the propagation of the crack is different for the three 
modes. 



Figure 4: Substructuring of the laminated composite structure 



2.1.1 Substructured formulation of the reference prob- 
lem 

The laminates structure E is decomposed into substructures 
and interfaces as represented in Figure ([3|. Each of these 
mechanical entities possesses its own kinematic and static 
unknown fields linked by its behavior. The substructuring 
is driven by the will to match domain decomposition inter- 
faces with material cohesive interfaces, so that each substruc- 
ture belongs to a unique ply P and has a constant linear 
behavior. A substructure E defined in Domain is con- 
nected to an adjacent substructure E' through an interface 
^EE' = dO^E^dO^E' (Figure Q). The surface entity Tee' 
applies force distributions F^, E_e' well as displacement 
distributions W^^, W_p^i to E and E' respectively. We write 

On a substructure E such that fl (5^1/ U (9^1^^) 7^ {0}, the 
boundary condition (U^jF^) is applied through a boundary 
interface F^^ . 

Let be the Cauchy stress tensor and £^ the symmetric 
part of the displacement gradient in substructure E. 

Then, the substructured quasi-static problem consists in 
finding s = {sE)EeE at a given step of the time integration 
scheme, where se = {ue->W.e->^e''—^^' which verifies the fol- 
lowing equations: 



• Kinematic admissibility of Substructure E: 
at each point of F^ , Uj, 



(9) 



Static admissibility of Substructure E: 

Tr (^^ = j^^ f^.UE*dQ. 

Jte 



(10) 



3 





Cohesive 
interfaces 



Perfect interfaces 



Figure 3: Substmcturing of the laminated composite structure 



• Linear orthotropic behavior of Substructure E\ 

at each point of Q^e^ ^-(-e) ^) 

• Behavior of the interfaces F^^/: 



at each point of T^e' ^ , 

^EEiWE^WE^^FE,FE.)=0 



(12) 



• Behavior of the interface at the boundary: 

at each point of F^^ , ^e^ {We ,Le)=^ (13) 

The formal relation ^ee' = ^^d ^Ej = named "inter- 
face behavior" can be made explicit in the two cases that we 
have to handle: 



• Perfect interface: 



• Cohesive interface: 



Fe^Fe^=0 
We-We^=0 



Le^Le'=^ 
s^pp'{Wei-We^Fe)=Q 



(14) 



(15) 



where Substructure E (respectively E') belongs to Ply P 
(respectively P'), such that P <P' . 

In the same manner, and for instance in the case of a pre- 
scribed displacement boudary interface Te^ , the formal rela- 
tion ^E^ = reads: 

We = U^ (16) 



2.2 Two-scale iterative resolution of the sub- 
structured problem 

2.2.1 Introduction of the macroscopic scale 

In order the scalability of the method to be ensured, a global 
coarse grid problem is solved at each iteration of the solver. 
The definition of the macroscopic fields required to construct 
this problem is done on the interface only. 

At each interface F^^/, the interface fields are split into a 
macro part ^ and a micro part the former belonging to a 
small-dimension subspace (e.g. 4 macro degrees of freedom 
per interface in 2D, 9 in 3D). 



Fe 



7M 



(17) 



Given the macrospaces and on Interface F^, the 
unicity of the decomposition of the interface fields in macro 
and micro data is ensured by uncoupling the interface virtual 
works: 

(18) 

Usually, a common macro basis for both the kinematic and 
static interface macro fields is chosen. Numerical tests 
showed that in order to ensure the numerical scalability of 
the method the macro basis should extract at least the linear 
part of the interface forces (see Figures ([5]) and (|7])). Indeed, 
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Figure 6: Substructuring of the laminated composite structure 

this macro space contains the part of the interface fields with 
the highest wavelength. Consequently, according to the Saint- 
Venant principle, the micro complement (which, as explained 
in next subsection, is found iteratively through the resolution 
of local problems) has only local influence. 

2.2.2 The iterative algorithm 

The iterative LaTIn algorithm, which enables the resolution 
of nonlinear problems, is here applied to the resolution of the 
substructured reference problem with nonlinearities localized 
in the interfaces. 

The equations of the problem are split into two groups: 

• linear equations in substructure and macroscopic inter- 
face variables: 

• static admissibility of the substructures 

• kinematic admissibility of the substructures 

• linear behavior of the substructures 

• equilibrium of the macro interface forces 

• local equations in interface variables: 

• interface behavior 

The interface solutions s = (s^j^^E = {W-E^E^E)EeE to the 
first set of equations belong to Space A^, while the interface 
solutions s = (sE)EeE = {W.E^£-E)EeE to the second set of 
equations belong to F. The converged interface solution Sref 
is such that: 

SrefeAdpr (19) 

The resolution scheme consists in seeking the interface so- 
lution Sref alternatively in these two spaces: first, one finds a 
solution s„ in A^, then a solution in F. In order for the 

two problems to be well-posed, search directions and E~ 
linking the solutions s and ? through the iterative process are 
introduced (see Figure ([6])). 

Hence, an iteration of the resolution scheme consists of two 
stages: 



• a local stage: 

Finds^^^ e r such that (^s^^+i —^n^ ^ (20) 

• a linear stage: 

FindSn^i G A^ such that {sn^i — s^^+i^ ^ E~ (21) 

In the following sections, the subscript n will be dropped. 

Local stage In the local stage, local problems are solved at 
each point of the interfaces T^^i : 

Find {F_^ , , F^/ , W_^f ) such that: 

^ee'{We,We.1e1e')=0 (22) 
{lE-FE)-k^{WE-WE)=0 
{Ie'-Fe')-4{We^-We.)=0 

the two last equations of this system being the search direction 
. In the case of a cohesive interface. Problem ([22]) is non- 
linear, and solved by a modified Newton-Raphson scheme. 

Local linear problems are also solved at each point of the 
boundary interfaces F^^ : 

Find {F_E , W_E)such that: 
r I^eA^e,Ie)=Q (23) 

I {lE-FE)-kUt.E-W.E)=^ 

Linear stage The linear stage consists in solving linear sys- 
tems on each substructure under the constraint of macroscopic 
equilibrium of the interface forces. 

on interface Tee' , F_e^E-e'=^ (24) 

Macroscopic admissibility of displacements could also be en- 
forced. In the case of perfect interfaces, it would be easy to 
derive. In the case of non-homogeneous or nonlinear behavior 
at the interfaces, the macro condition would not be practically 
(and sometimes theoretically) feasible. 

Condition ([24]) is incompatible with the monoscale search 
direction E~ coupling the interface displacement and forces 
fields at the linear stage, which reads: 

on interface F^ , {Le-Le)^ {We-W_e)=0 (25) 

Hence this search direction is weakened and verified at best 
under the macroscopic constraint. Technically this is realized 
using a Lagrangian whose stationarity leads to a modified lo- 
cal search direction: 

V WE^'eWE, I {Fe-Ie).We'' dT 

+y {kE{WE-WE)-kEW ywE^'dr^o 
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Figure 5: The linear macro basis for a plane interface 



and to an equation expressing the continuity of the macro- 
forces: 



E ^^E 



where W_ is unique for Interface T^^i and set to zero on duO.. 
A way to solve this set of equations consists in introducing 

a relation coupling and W_ into ([27]). This relation is de- 
rived from the local equilibrium of each substructure ([T0| and 
from the local modified search direction ([26]). The problem to 
be solved for Substructure E becomes: 

/ TT{£{uE)K£{u/))d^^ [ ^^W^.W/jr 
JQe ^ -^^e ^^^^ 



where F_e = £-e~^ ^e ^e • 



Equation ( [28] ) is linear. Therefore, one can write a linear 
relation between the interface displacements and the loading: 

Jte 



(29) 



Operator is the dual Schur complement of Substructure E 
modified by the search direction, while Wk^ results from the 
condensation of the volumic loading on interface F^. 

The corresponding interface forces are obtained through 
the modified search direction ( [26] ) and projected onto the 
macro space: 



I F^.w^'^dT^ I (l: 



where 



JTf 



Ee W. 



^ {FE-kE (meFe^w^e^)) .w'^'^dr 



]Le is classically viewed as a homogenized behavior of Sub- 
structure E and is calculated explicitly for each substructure 
by solving local subproblems ( [28] ) taking the vectors of the 
macro basis as boundary conditions on F^. 

This relation is finally introduced into the equation express- 
ing the admissibility of the macroforces ( [27] ), leading to the 
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so-called macro-problem: 



E -^^E 



M -Mi. 



dT- 



eJ^e 



dT 



(31) 

The macro-problem is discrete by nature. Hence, its an al- 
gebraic form = F^, where is the vector of the 
components of the Lagrange multiplier in the macro basis, is 
also used in the following. 

The right-hand side of Equation ( [3T] ) can be interpreted as 
a macroscopic static residual obtained from the calculation of 
a single-scale linear stage. In order to derive this term, the 
problem ([28]) must be solved independently on each substruc- 

-M 

ture, setting W_ to zero. The resolution of the macroscopic 
problem |3T] ) leads to the global knowledge of Lagrange mul- 
tiplier W_ , which is finally used as prescribed displacement 
to solve the substructure independent problems ( [28] ). 

In order to perform the resolutions of ( [28] ) on the substruc- 
tures, finite element method is used. Since the behavior of the 
substructures is linear, the stiffness operator of each substruc- 
ture can be factorized once at the beginning of the calculation 
and reused without updating throughout the analysis, which 
gives the method high numerical performance. 

Algorithm p.l| ) sums up the iterative procedure which has 
been described in this section. 



Algorithm 2.1 The 2-scale domain decomposition solver 
1: Substructures' operators construction 
2: Computation of the macroscopic homogenized behavior 

on each substructure 
3: Global assembly of the macroscopic operator 
4: Initialization Sq G F 
5: for n = 0,...,N do 
6: Linear stage: computation of s„ G A^ 

□ Computation of the macroscopic right-hand term 

M 
E 

on each substructure 

□ Global assembly of the macroscopic right-hand 
term 

□ Macro problem resolution 

□ Micro problem resolution 
Local stage: computation of s^^ i 

□ Boundary interfaces Te^ 

□ Internal interfaces T^e' 
Calculation of an error indicator 

end for 




Propagation of the delamination crack 



Initially-delaminated interface 



Figure 7: Four-ply DCB test case 



2.3 Delamination analysis example 

A first example of quasi-static delamination analysis is shown 
in Figure ([7]). The problem consists in a [0/90]^ double can- 
tilever beam (DCB) case. The loading leading to mode I 
quasi-static crack's propagation is increased linearly over 10 
time steps. The first three of them correspond to the initiation 
of the delamination and the remainder to the crack's propaga- 
tion. 

This assessment is realized with a C++ implementation of 
the mixed domain decomposition method capable of handling 
the quasi-static analysis of 3D nonlinear problems. The paral- 
lel computations use the MPI library to exchange data among 
several processors. 

Each processor is assigned to a set of connected substruc- 
tures (along with their interfaces); it calculates the associated 
operators and solves the local problems. This tends to reduce 
the number of interfaces duplicated among several processors 



(Figure (17)) and is achieved technically through a METIS 
routine. 

This very simple test case already exhibits numerical diffi- 
culties: 

• The convergence rate of the LaTIn-based strategy is 
highly dependant on the search direction parameters. In 
the case of cohesive interfaces the iterative solver can 
even stagnate when using too small values for the search 
direction parameters. Hence, we describe in next section 
the way we set them in order to ensure convergence. 

• The method loses its numerical scalability when the 
crack's tip propagates. This phenomenon appears clearly 



in Figure ( [15] ) ("No sub resolution" label). When the de- 
lamination process propagates (Time steps 4 to 10), the 
number of LaTIn iterations to convergence becomes very 
large. A solution to this problem is developed in Section 
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3 Analysis of the iterative algorithm be solved directly, 
parameters 

The search direction parameters {k'^)(EeE) (^^)(£gE) 
introduced as positive definite symmetric operators. It has 
been empirically shown in previous studies that an optimal 
set of these operators exists. Though, these optimal operators 
are known to be difficult to interpret theoretically, especially 
when using complex interface behaviors. In addition, they 
can be expensive to compute even in simplified cases (perfect 
interface behavior). Our goal is here to find an efficient local 
approximation of the search direction operators for debonding 
analysis. 




3.1 Search direction 

Using too small a value for Parameters {k'^)(EeE) can lead to 
the stagnation or divergence of the algorithm. Actually, in 
this case, the interface solution is seeked at the local stage 
in a truncated space F, the solutions in the softening part of 
the local cohesive behaviour being unreachable. Figure ([8| 
illustrates schematically this idea. The constitutive law of a 
cohesive interface and the search direction equation are 
projected on the normal direction to the interface, at a given 
integration point. The solution seeked at the local stage is 
the intersection of the two curves obtained. This solution is 
computed by a Modified Newton-Raphson scheme. It clearly 
appears here that the part of the constitutive law within the 
dotted frame cannot be reached when using this iterative pro- 
cedure. Thus, the global minimum of the problem may not be 
found through the LaTIn iterations. 




Figure 8: Using too small a value for Search direction param- 
eter 

We thus choose to set {k'^)(EeE) to infinite values (see Fig- 
ure [9]) and focus on the choice of Search direction E~ only. 
Although slightly improved convergence rate could be ob- 
tained using classical conjugate search directions and E~, 
the time costs of the local stage drop as the local problems can 



Figure 9: Using infinite value for Search direction parameter 



3.2 Search direction E : interpretation 

The search direction E~ must be separated into a macro part 
E~^ and a micro part E~™ in order to be interpreted. Equa- 
tion ([26]) can be rewritten: 

* e#f, ^ [{FE-lE)).wfdr 



w^-w^] .w¥* dr = o 



(32) 

G^-i", / ({FE-FE)).wTdr 

r / ^ X (33) 

+ J^^[k-E'" {We -We)) .WT dT = Q 

where is the space orthogonal to with respect to the 
L^iTE) inner product. 

Perfect interfaces Previous studies flTl [211 have focused 
on the choice of the micro search direction parameter 
{kE^)EeE- Classically, when dealing with perfect interfaces, 
the optimal parameter k^"^ on interface Te can be linked to 
the Schur complement of the structure occupying the domain 
Q,\Q,E- As the introduced microscopic interface quantities 
have a local influence, this Schur complement can be calcu- 
lated in the vicinity of Interface Te . In the case of isotropic 
and homogeneous materials. The scalar £^/L has been shown 
to be a good approximation of this local operator (211, E be- 
ing the Young's modulus of the adjacent substructure and L 
a characteristic length of the interface. As explained in para- 
graphs, computing this microscopic optimal search direction 
parameter is not of primal interest in our case. 

The {k^^)EeE search direction parameter can be inter- 
preted by considering an interface Tee' separating two ad- 
jacent substructures E and E\ In order to simplify the expla- 
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nations, we consider here a unique search direction parame- 
ters on Tee' for both substructures. Taking into account the 
macroscopic equiUbrium of the interface forces along the in- 
terface r^£/, on can derive from ([32]) : 




This equation means that k^^, is a macroscopic stiffness of 

_M 

interface Tee'- In the case of perfect interface, setting k^^, 
to the infinity equals to enforcing the macroscopic displace- 
ment continuity through the interface. As the the macroscopic 
equilibrium of the interface forces is ensured by equation ( [27| ), 

_M 

this choice of the search direction parameter k^^, leads to the 
complete enforcement of the interface macroscopic behavior 
at the linear stage : 

f ywfewi^ I {FE^FEf).wf dr = o 

I r^E' (35) 

vw^'^G^/ / {WE-WEf).wf dr = o 

Though, as the substructures are small compared to the 
size of the structure, the characteristic length introduced 
to approximate the microscopic optimal search direction 
parameter is very small. Hence, setting k~^, to is 
close to ensuring the macroscopic displacement continuity 
through perfect interfaces. Moreover, as it will be explained 
further on, the choice of a unique search direction paramater 
increases the simplicity and efficienty of the iterative strategy. 

In the case of more complex interface behaviors, like soft- 
ening, ignoring the influence of the macroscopic stiffness 
(^^^)(£eE) can lead to non physical solutions or even to the 
divergence of the iterative process. This idea is illlustrated by 
the results in Figure ([TO]). The test case presented is a four- 
plies composite beam, the black parts of the cohesive inter- 
faces corresponding to initial delamination. The third inter- 
face is artificially weakened. The solution labelled "Refer- 
ence solution" is obtained by a modified Newton algorithm, 
while the two following ones result from a LaTIn computa- 
tion, two different set of search direction parameters being 
used. 

Thus, an effort must be made to correctly set the interface 
conditions prescribed at the linear stage. In the case of simple 
interface behaviors (including the perfect interface behavior 
studied in the previous paragraph), a "physical" macroscopic 
interface condition can be prescribed beetween sub-structures 
at the linear stage. We will derive the general case from the 
analysis of these specific studies. 




Figure 10: Converged solutions reached by the LaTIn algo- 
rithm using different sets of search direction parameters 

Homogeneous isotropic elastic interface Let us set the 

_M 

macroscopic search direction parameter k^^, of an elastic in- 
terface r^^/ to 2^^(1 — J), where is the local scalar stiff- 
ness of the interface. Equation ([34]) now reads : 

{V W^GyT/ / {Fe^Fe,).W_T dT = 
/ {FE-k\\-d){WE.-WE)).wf dT = Q 

(36) 

Hence, setting ^k^^, to the interface stiffness equals to pre- 
scribing the macroscopic interface behavior as an interface 
condition at the linear stage. 

Delaminated interfaces under traction In the case of a 
completely damaged cohesive interface loaded with traction 
(so that the gap of interface displacements is positive), the 
converged interface force fields are null. Thus, both micro 
and macro search direction parameters should optimaly be set 
to zero. Indeed, F^e — E_e' = on interface Tee' these 
quantities result from the local stage and verify the interface 
behavior. Consequently, setting k~^, = = leads to the 
enforcement of the relation = F^/ = (equations ([32]) and 
([33])). 

General case We can conclude from these three simplified 
cases that, as the macroscopic equilibrium of the interface 
forces is enforced at the linear stage through the choice of 
Admissibility space Ad, one can ensure the verification of the 
complete linearized interface behavior by the macroscopic in- 
terface fields at this stage by setting the E~ interface param- 



9 




1 23456789 10 

Time step 

Figure 1 1 : Use of small search direction parameters for delaminated interfaces 
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eters to specific values. Thus, together with the choice of a 
rigid E+ search direction, this choice leads to a hybrid itera- 
tive strategy : 

• a modified Newton-Raphson scheme on the macro-part 
of the solution 

• a Latin-type algorithm on the micro part, the search di- 
rections and E~ being non-conjugated (which is not 
classical) 

Though, when dealing with cohesive interfaces with non- 
constant stiffness, with interfaces partially delaminated, or 
with delaminated interfaces loaded with both local traction 
and local compression, the interface behavior that should be 
verified by the macroscopic fields is not explicit. Yet a too 
coarse approximation of this behavior might drive the algo- 
rithm towards a local minimum of potential energy instead of 
a "more physical" global minimum. 

Moreover, the mechanical interpretation of the search di- 
rection parameters {k^)(EeE) requires to introduce a micro 
search direction E~™ and a macro search direction E~^. But 
when using this search direction separation, we found out 
that the CPU time increased, because of the large amount of 
projections of the interface fields in the macrospace required 
through the iterative process. 

A practical way to choose a common micro and macro 
search direction parameter that ensures a good interface con- 
dition is the subject of next subsection. 

3.3 Search direction E~: practical choice 

The {k^)(EeE) search direction parameters are set with respect 
to the interface behavior, as explained bellow: 

• perfect interfaces: {k^)(EeE) set to the optimal micro 
values described previsouly. As the characteristic length 
L of the interface is very small, these values are high 
enough to ensure the continuity of the macrodisplace- 
ment through the perfect interfaces. 

• undamaged or partially damaged cohesive interfaces: 
{k^)(EeE) set to the optimal macro values for initially 
undamaged interfaces. In the case of an interface F^^/, 
we thus choose: 

/ 2^0 \ 

k-j,, = 2^0 (37) 

V A...,.) 



branching to non-physical solutions has been observed 
in our cases. 

• delaminated interfaces: {k^)(EeE) set to the same 
values used for undamaged interfaces, unless the contact 
status of the interface is known (i.e.: unless all integra- 
tion points are in compression, or unless all integration 
points are in traction). In the last cases, the search direc- 
tion parameters are updated as follows : 

• delaminated interface under compression 

/2^o 0\ 
kEE'=\ (38) 

\ / (nEllll) 

• delaminated interface under traction 

/:-^,=OxId (39) 

Obviously, updating E+ with respect to the contact sta- 
tus of the cohesive interface requires a re-assembling 
step of the macroscopic global operator. Potentially, this 
method can be expensive, unless the macroscopic prob- 
lem is solved using a parallel strategy (see Section |5]). 

The results of this procedure for the cubic test case rep- 
resented Figure ^ are shown in Figure ([TT]). All the 
interfaces between adjacent sub- structures are granted a 
cohesive behavior. The prescribed loading leads to initi- 
ation of the delamination (time steps 1 to 2), opening of 
the cracks (time steps 3 to 7) and closing of the cracks 
(time steps 8 to 10). In the first case, the search direction 
parameters {k'^)(E^E) constant throughout the anal- 
ysis. In the second case, they are updated as explained 
previously, the contact status of the delaminated inter- 
faces being checked every ten LaTIn iterations. Clearly, 
the number of iterations to convergence is significantly 
reduced when delamination occurs, the delaminated in- 
terface being under traction. 

• prescribed forces ( respectively displacement) interfaces: 
{k^)(EeE) ^re set to a very low (respectively high) value 
so as to enforce the boundary condition through penal- 
ization to the adjacent substructure. 



even though it cannot be shown theoritically, this strat- In the end, it comes out that the retained parameters are 

egy has always led to the convergence of the iterative al- quite independant on the shape of the cohesive law. Thus, 

gorithm. Moreover, as the interface conditions between good performance results can be expected from other damag- 

substructures are connected to mechanical properties, no ing interface models, without significant adaptation. 
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Figure 13: Extraction of a subproblem 
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4 Subresolutions in the crack's tip 
vicinity 

The drop in the convergence rate occurring when the cracks 
propagate can be explained by two main phenomena: 

• the singularity near the tip of the crack is very poorly 
represented by linear macro quantities (see Figure ([12])). 
Therefore the complementary "micro" parts of these 
phenomena, which are calculated iteratively through the 
resolution of local problems, have an influence on a sig- 
nificant part of the structure. In order to maintain the 
scalability of the method, a method to filter this global in- 
fluence from the "micro" quantities in the process zone 
and then transmit it to the whole structure must be de- 
signed. 

• the prediction of the location of the crack's tip at a given 
time step requires the quasi-convergence of a large num- 
ber of consecutive equilibrium states. Consequently, the 
propagation of the crack is very slow as the iterations 
proceed. 

A first solution would be to enlarge the macro space so that 
the totality of the numerical influence of the crack would be 
systematically transmitted over the whole structure. Because 
such a strategy would imply to update (reassembly and refac- 
torization) the macro problem at each evolution of the crack, 
it is not computationally realistic. However, the very-large- 
variation-length part of the solution is determined correctly in 
most of the structure since the very first iterations. Therefore 
as exposed in this section, one can choose to solve "exactly" 
the highly nonlinear problem in the crack's front vincinity at 
each global LaTIn iteration. 

4.1 Principle of the sub-resolution strategy 

The technique consists in extracting a part Q^sub of Domain Q,. 
The converged solution of the extracted nonlinear subproblem 
on Qsub is sought using the two-scale domain decomposition 
strategy described in Section [2] (see Figure ([13])) along with 
Algorithm (4.1 )). This idea is similar to the concept of non- 
linear relocalization developed in [5) [271, in which the authors 
used a domain decomposition method and performed a non- 
linear analysis in each substructure after a resolution of the 
condensed global linearized problem. In our case, the non- 
linear relocalization is carried out on a set of substructures 
because the nonlinearities are localized at the interfaces of the 
domain decomposition scheme. 

Let us concentrate on two main difficulties: 

• the choice of the boundary conditions for the subproblem 



• the choice of the size and position of the extracted sub- 
domain 

4.2 The boundary conditions of the subprob- 
lem 

The subresolution is carried out at each step of the global it- 
erative resolution, which results in unbalanced substructure 
and interface fields. Numerical tests showed that the pre- 
scribed conditions applied on the boundary dQ^sub of domain 
Qsub must be in Space (i.e. they must result from a lin- 
ear stage of the resolution). Indeed these fields are in global 
static equilibrium over the whole structure, whereas the so- 
lutions in Space F are in equilibrium only locally and do not 
self-equilibriate O^sub- 

Thus, subresolutions can be interpreted as an enhancement 
of the linear stage. In order the process zone to keep matching 
the remaining of the structure, Robin boundary conditions are 
prescribed on dQsub using search direction E~ as interface 
stiffness parameter. 

4.3 Adaptivity of the subproblem 

In order to extract the subproblem automatically, we choose to 
select a set of substructures and interfaces in a box surround- 
ing the interface with the highest damage rate at the end of the 
global step (see Figure ([13])). 

The influence of the size of the extracted subdomain is 
shown in Figure ([14]). 

4.4 Results 

The resolution of a subproblem around the crack's tip leads to 
a convergence rate of the global resolution which is indepen- 
dent of the time step of the analysis (i.e. independent of the 
area of the interface which becomes delaminated in one time 
step), which means that the numerical scalability is restored. 
As a result, the local inversion time for the problem shown 



in Figure ( 13 ) (with the smaller subiteration domain) was cut 
in half. This estimate does not take into account the fact that 
the macroproblem is much smaller in the case of the subitera- 
tions; furthermore the gain increases as the ratio between the 
size of the process zone and the size of the structure decreases. 

Thus, this method can lead to a reduction in the number 
of calculations. However, this reduction would be ineffective 
unless the subproblem is re-parallelized. Indeed, using the 
initial allocation among the parallel processors would adress 
the extracted subproblem to only a very small number of pro- 
cessors. Another solution, easier to implement but potentially 
less efficient, would be to perform independent subiterations 
systematically on all the processors. 
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Algorithm 4.1 The subresolution strategy algorithm 
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Figure 14: Influence of the size of the subproblem 
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Figure 15: Subiterations around the crack's tip 
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Substructures' operators construction 

Computation of the macroscopic homogenized behavior 

on each substructure 
Global assembly of the macroscopic operator 
Initialization Sq G F 
for n = 0,...,N do 

Linear stage: computation of G 

Local stage: computation of i G T 

Subproblem extraction 

□ Location of the substructures requiring subresolu- 
tion 

□ Application of mixed boundary conditions 

□ Assembly of the macro subproblem 
for J = 0, . . . , m do 

Subproblem linear stage 
if j <m—l then 

Subproblem local stage 
end if 

Local error indicator 
end for 

Local stage on the boundary interfaces of the subprob- 
lem 

Calculation of an error indicator 
end for 



5 Third-scale resolution 

The substructuring described in Section |2] results in a very 
large macro problem and in an unnecessarily refined macro- 
scopic solution. In order to solve large problems such as the 



one represented Figure ( 16 ), we need to focus on: 

• the parallel resolution of the macroproblem, 

• the selection and transmission of the large-wavelength 
part of the macro solution. 

These two elements can be introduced into the method 
through the use of any Schur-complement-based domain de- 
composition method 1 141 . We chose to implement the BDD 
method 1221 to solve the macroproblem. 

5.1 The balancing domain decomposition 
method for the macroproblem 

5.1.1 Partitioning of the macroproblem 

The substructures of the initial partitioned problem are 
grouped into super- sub structures E separated by super- 



interfaces r^^/ (Figure (17)). Practically, each super-sub- 
structure is made of the whole set of sub- structures assigned to 
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a given processor of the parallel computing architecture. The 
algebraic problem to be solved within each of these super- 
substructures reads: 



jMiE) 
^ bi 



tm(^) 
^ bb 



-If) 



(40) 

where the ^ superscript has been omitted, the subscripts b 
and / refer to the super-interface quantities and to the inter- 
nal quantities of the super-substructures respectively. B^^^ 
and A^^) are signed Boolean localization operators. The sec- 
ond equation of System ( |4Q| ) expresses the continuity of the 
kinematic unknowns, while the third equation expresses the 
equilibrium of the nodal reactions at the interfaces between 
super- substructures. 

The local equilibria are condensed at the super-interfaces 
through the introduction of Schur complements S^^^ and con- 



densed forces F 



(E). 



(E) 



AE) 



(41) 



where 



jMiE) 

bb ' 



AE) AE) _ 



T M(^) T M(^) 

bi ^ a 

T M(^) t M(^)~ 



jUiE) 
^ ib 

'Ae) 



The continuity of displacement is achieved automatically 
through the introduction of a unique super-interface macrodis- 
placement W^. Then, the continuity equation of the interface 
reaction forces yields: 



where < 



E 



(42) 



5.1.2 Resolution of the super-interface problem 

The condensed macroproblem is solved iteratively through a 
conjugate gradient algorithm. Classically, this resolution re- 
quires only matrix- vector products and dot products, which 
are compatible with parallel architectures. 

The recommended preconditioner for a parallel use of this 
algorithm is what is called the Neumann preconditioner: 



(43) 



Algorithm 5.1 Projected preconditioned conjugate gradient 



Initialize W^, = {VWt J ^ C{C^ SC)'^ 
Calculate ro = iv — SW^q 
Calculate zo = PS~^ro and set wq = zo 
for 7 = 0, . . . ,m do 



a, 



{rj,Zj)/{Swj,Wj) 



end for 



-{Swj,Zj^i)/{wj,Swj) 



S*^^) being a generalized inverse of the Schur complement of 
Super- substructure E. 

The use of this preconditioner means that the inverse of the 
global super-macro operator is approximated by the assembly 
of the inverses of the local Schur complements. Let us note 
that the description chosen for the interface macrofields pre- 
cludes the existence of degrees of freedom belonging to more 
than two substructures; consequently, no scaling is required 
in the preconditioner (at least when the interfaces are not too 
heterogeneous). _ 

Since the product S~Vj+i consists in solving Neumann 
problems for each super- sub structure E under the loading 
r^+i , one must ensure that r^+i is self-balanced in the sense of 
E. Therefore, we introduce a projector P which projects the 
residual onto the space orthogonal to the kernel of the super- 
substructure at each iteration of the conjugate gradient. 

Thus, the solution is sought in the form: 



(44) 



where 

C^SP 



%o =C(C^SC)-^ C^Fc 
P = I-C(C^SC)^C^S 



Matrix (C^SC) corresponds to a coarse representation of 
the global stiffness of the structure. Operator C must contain 
at least the rigid body modes of the super- sub structure. Then, 
the initialization T^^q is achieved as a combination of the rigid 
body modes of the local stiffness operators (range(C)), while 
the remaining part is sought iteratively in the supplementary 
subspace ker(C^S) through the projector P. 



5.2 Results 



Figure ([18]) shows the convergence rate of the LaTIn algo- 
rithm when the conjugate gradient scheme for the condensed 
macroproblem is stopped after a fixed number of iterations. 
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Figure 18: LaTIn convergence curves for different numbers 
of macro iterations 



The test case is the holed plate under traction loading rep- 
resented in Figu re (p^ with the super- substructuring pattern 
given in Figure (|17|. The structure is divided into 520 sub- 
structures, separated by 1350 interfaces. The number of mi- 
cro (respectively macro) degrees of freedom involved in this 
problem is 3.4 x 10^ (respectively 12150). It appears clearly 
that only very few iterations of the conjugate gradient scheme 
are required to get the necessary part of the macrodisplace- 
ment Lagrange multiplier leading to the multiscale conver- 
gence rate of the LaTIn algorithm. Hence, high accuracy of 
the macroscopic resolution is not necessary to transmit perti- 
nent piece of information on the whole structure. Typically, 
the algorithm is stopped when the residual error (normalized 
by the initial error) falls below 10~^. The admissibility of the 
macroforces is thus enforced on a third level, which is suffi- 
cient to determine the part of the solution which needs to be 
transmitted at each iteration of the resolution. 

Figures ([19]) and (20) show how well the method scales 
both in computational time and memory usage when em- 
ployed on modern hardware parallel architectures (distributed 
memory clusters). 



6 Control of the loading sequence 

The incremental version of the LaTIn algorithm, like Newton- 
Raphson algorithm encounters numerical difficulties when 
used to carry on analysis beyond global limit-points or snap- 
backs. 

This section focuses on the discretized assembled nonlinear 
problem : 

m<n 

)Un=F„ (45) 
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Figure 19: LaTIn convergence curves for different numbers 
of macro iterations 



Global equilibrium states are sought successively at each time 
step n {0 <n<N) using an implicit time integration scheme. 
The nonlinear equilibrium problems coupled with a local arc- 
length control are solved by a modified Newton algorithm, the 
linear prediction steps being handled by the 3-scale domain 
decomposition strategy, as described in next subsection. 

The control algorithm is activated locally during the time 
analysis when the LaTIn strategy fails to converge. Con- 
versely when the control algorithm results in a succession 
of re-increasing loadings, the solution algorithm is switched 
back to the LaTIn method. 



6.1 Local control 

As usual in arc-length methods |[28|[6]|, the amplitude of the 
loading A„ is linked to the global displacement in such a way 
that the norm of the {AUn^AXnF) takes a predefined fixed 
value, the A . unknowns being the increments of these quanti- 
ties between time steps n — l and n. In our case, these global 
unknowns are not of primary interest (301 [H. Instead, let us 
introduce the control equation as: 



c{Un)AUn=Al 



(46) 



where c{U) is a Boolean operator extracting the maximum 
value of the local displacement gap over all the cohesive inter- 
faces of the model, and A/ is a given value. Thus, the loading 
increment A is controlled by a local variable which is closely 
related to the maximal local damage increment of the struc- 
ture. 
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^ c{Ul)K{Ul)-^l 



(50) 



Then, Operators K{Un) and c{Un) are updated with respect 
to the kinematic field Ui^i found at the prediction stage of 
the modified Newton algorithm ( [50| ) and, unless the residu- 
als of the updated equilibrium and control equation are small 
enough, a new iteration is performed. The norm of the resid- 
ual of the control equation is not used to stop the Newton iter- 
ations. In facts, the given value of the maximum local damage 
increment has no physical meaning. Yet, using this numeri- 
cal technique ensures that the evolution of the loading permits 
to follow the global behavior of the structure. Consequently, 
any converged equilibrium state found can be used to perform 
a new time-step computation. 



6.2 The arc-length resolution 



6.3 Parallel calculation 



Thus, at each time step n, we seek the solution (Un^^n) of the 
following discrete system: 



f{Un,Xn)=K{Un)Un-XnF = 
g{Un,Xn)=c{Un)AUn-M = 



(47) 



The first-order expansion of the equilibrium equation around 
the point {Ui^Xi) yields: 

dU HUM 
+ ^ , ^ (Ai+i-Ai)=0 

dx 

(48) 

We use a modified Newton algorithm, which means that the 
tangent operator |^|(^ is approximated by K(t/^), while 

W\{u' A') approximated by c(t/^). Using the relations ex- 
pressing that System ( |47| ) is verified at time step (^ — 1), one 
gets: 

1 c{U^)AUi-^'=Al 



(49) 



The introduction of the linearized equilibrium into the lin- 
earized control equation leads to the expression of the loading 
parameter at Iteration (/+!), then to the displacement solu- 



Our attempt to solve the linear problem ( [47] ) has been unsuc- 
cessful. This can be explained by the fact that the control 
equation is global over the whole structure and non-linear. 
Thus the classical separation of the linear equations on one 
hand and the local non-linear equations on the other hand, 
which is the basic idea of the LaTIn method, cannot be made. 
Nevertheless, the prediction step of the Newton algorithm de- 
scribed in the previous section requires the resolution of the 
linear system K{U^)~^F. We propose to use the LaTIn mixed 
three- scale domain decomposition method to find the solution 
to this linear system. The method is still efficient for two rea- 
sons : 

• The non-linearities being computed are still localized 
on the interfaces of the domain decomposition method. 
Consequently, no re-assembling step is required as the 
Newton iterations proceed. 



• using a full nonlinear LaTIn solver or using a Riks solver 
with a parallel LaTIn resolution of the prediction step 
requires very similar computations, which means that 
switching from the LaTIn algorithm to the arc-length al- 
gorithm is straightforward 

Nevertheless, one should notice that the linearity of the pre- 
diction step of the Riks solver makes the subiteration tech- 
nique described in Section ^ non-relevant (or at best less 
effective). 
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Figure 23: Behavior of the holed plate under traction 
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Figure 21: Influence of the initialization of the LaTIn linear 
solver on the number of LaTIn iterations required to achieve 
convergence at each Newton iteration 
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Figure 22: Influence of the stopping criteria on the numbers 
of iterations of the Newton algorithm and of the LaTIn linear 
solver, to compute one Newton increment 



6.4 Numerical improvement of the Newton al- 
gorithm using the LaTIn method as a lin- 
ear solver 

Two elements can reduce the cost of the arc-length algorithm 
significantly: 

• the initialization of each iteration of the local arc-control 
algorithm with the interface quantities found after con- 
vergence of the LaTIn resolution scheme in the previous 
iteration. The number of LaTIn iterations decreases very 
rapidly as the iterations of the Newton algorithm pro- 
ceed because the changes in the secant operators c{Un) 
and K{Un) become less and less significant. On Figure 
( [2T] ), the two first time steps of the holed plate test case 
(see Figure ([16])), which correspond to different levels of 
non-linearity, are computed using a stopping criteria of 
the LaTIn linear solver set to a given value. The drop 
of the number of LaTIn iterations required as the New- 
ton algorithm goes on appears clearly in these two cases. 
As this idea systematically improves convergence for no 
extra cost, it is now used by default. 

• the crossed-optimization of the stopping criteria of (non- 
linear) Newton solver and (linear) LaTIn solver. Basi- 
cally, the LaTIn algorithm could be converged to a very 
low value of the residual of the linear system at each pre- 
diction step of the Newton scheme. Though, our tests 
show that the first iterations of Newton, leading to a high 
value of the residual of the non-linear system, do not re- 
quire an exact resolution of the prediction step. In order 
to illustrate this idea, the number of Newton iterations 
for one time step of a DCB test (Figure ([TSjtop)), and the 
associated total number of (linear) LaTIn iterations are 



plotted in Figure ( 22 ) as functions of the ratio of these 
two errors . This shows clearly that in order to use the 
method most efficiently (i.e. with the smallest total num- 
ber of LaTIn iterations), the convergence threshold of the 
LaTIn method should be set to an error very close to the 
current Newton error. 



6.5 Results 



Figure ( 23 ) shows the global force v^. displacement curve ob- 
tained for the holed plate test case ([16]) using the arc-length 
algorithm described above. The damage in the interfaces 
loaded in shear mode (interfaces [0/90]) is also represented 
at four equilibirum states of the time analysis. Several very 
sharp snap-backs appear in the global behavior curve of the 
structure, and are efficiently handled by this locally controlled 
Riks' algorithm. 
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7 Conclusion 

The accurate prediction of delamination in large process 
zones of laminate composite structures requires refined mod- 
els of the material behavior. Such descriptions lead to the 
resolution of huge systems of equations. In order to calcu- 
late the exact solution of such a refined model, we used a 
two-scale domain decomposition strategy based on an itera- 
tive resolution algorithm. This method is particularly well- 
suited for laminate models in which 3D and 2D entities are 
introduced separately. 

This strategy has been improved in order to make it capable 
of handling very large delamination problems. A systematic 
analysis of the features of the method at the different scales 
has been conducted. It has first been shown that the classical 
scale separation was insufficient in the high gradient zones to 
provide numerical scalability. We thus developed a subres- 
olution procedure which preserved the numerical scalability 
of the crack propagation calculation. This analysis has also 
proved that a third scale was required. The second-scale prob- 
lem is then solved using a parallel iterative algorithm, which 
enabled the fast transmission of the very-large-wavelength 
part of the solution. 

In order to perform the quasi- static analysis beyond the 
global limit-points resulting from the local softening behavior 
of the structure, we used an arc-length-type algorithm to con- 
trol the magnitude of the loading. We showed that the com- 
putation steps required when using this algorithm were very 
similar to those of the LaTIn technique. Therefore, switching 
from one algorithm to the other was very easy. 

In future developments, this 3D process zone analysis tech- 
nique should be associated with a plate model analysis, which 
would be sufficient to obtain the solution in the low-gradient 
zones. 
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